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Abstract. We present an N-body computer code - aimed at studies of galactic dynamics - with a CPU-efficient 
algorithm for a continuous (i.e. time-dependent) stellar mass-loss. First, we summarize available data on stellar 
mass-loss and derive the long-term (20 Gyr) dependence of mass-loss rate of a coeval stellar population. We then 
implement it, through a simple parametric form, into a particle-mesh code with stellar and gaseous particles. 
We perform several tests of the algorithm reliability and show an illustrative application; a 2D simulation of a 
disk galaxy, starting as purely stellar but evolving as two-component due to gradual mass-loss from initial stars 
and due to star formation. In a subsequent paper we will use the code to study what changes are induced in 
galactic disks by the continuous gas recycling compared to the instantaneous recycling approximation, especially 
the changes in star formation rate and radial inflow of matter. 
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1. Introduction 

During the last decade, many N-body simulations (e.g. 
Friedli & Benz 1993; Junqueira & Combes 1996) as well 
as theoretical works on two-fluid gravitational instabili- 
ties (Jog 1992, 1996) have underlined the importance of 
the gas mass fraction for the dynamics and evolution of 
galaxies. 

Gas crucially conditions star formation and evolution 
of large-scale instabilities (e.g. bars, spiral arms). In turn, 
star formation and large-scale flows induced by these in- 
stabilities influence mass and chemical profiles of galaxies 
(e.g. Martinet & Friedh 1997; Portinari & Chiosi 2000). 

The gas content in a given area is determined by the 
competing processes of star formation and stellar mass- 
loss and by spatial flows of gas. A coherent picture of 
galactic evolution must therefore consistently couple stel- 
lar and gaseous dynamics with star formation and stellar 
evolution. 

Star formation is, in large-scale models of galaxies, typ- 
ically implemented through simplistic parametrizations 
reflecting mainly the energy supply from young stars and 
the gas mass fraction locked in newly born stars. On the 
other hand, stellar mass-loss has received much less atten- 
tion from galactic N-body modellers since it was for a long 
time considered as a secondary issue. But nowadays, both 
observations and stellar evolutionary models (see Sect. 2 
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for references) indicate that the gas mass fraction resti- 
tuted by stars is huge and may reach, when integrated 
over the stellar mass spectrum, some 45% over the Hubble 
time. Stars thus represent not only a place for permanent 
gas blocking but also an important temporary reservoir of 
gas that will be gradually reinjected into the interstellar 
medium. 

Recent computer models of galactic dynamics usually 
approximate the stellar mass-loss as instantaneous, i.e. 
happening at the moment of stellar birth. Nonetheless, 
stars lose matter during all their lives. Stellar lifetimes 
span a very large range and in the case of low-mass stars 
they compare with or overpass the Hubble time. The use 
of the instantaneous recycling approximation, which was 
proposed by Tinsley (1980) for high mass-stars, is there- 
fore not satisfactory for the whole stellar mass spectrum. 

From the above emerges an obvious motivation for 
building a computer code able to follow the galactic struc- 
ture and dynamics together with a non-instantaneous gas 
recycling. One can expect that the connection of such a 
recycling with the dynamics will affect large-scale gravita- 
tional instabilities, spatial flows of matter, star formation 
and gas consumption rates, etc., and thus the long-term 
evolution of galaxies. 

The development and presentation of such a code is 
the central aim of this paper. Our model is innovative es- 
pecially in introducing a continuous (i.e. time-dependent) 
gas recycling scheme, grafted on an underlying N-body 
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code using standard techniques for computing gravita- 
tional field (particle-mesh method) and gas dynamics 
(sticky-particles) . 

The text is organized as follows. In Sect. 2, we de- 
scribe data and steps in deriving the time-dependence of 
mass-loss of a computer stellar particle. This curve, simply 
parametrized, is the key input information for our N-body 
code, however it can be useful also for other models (not 
necessarily N-body) of galaxies or star clusters. Section 3 
presents our computer code, the emphasis being put on 
the implementation of the mass-loss. An instructive sim- 
ulation of a barred galaxy is shown in Sect. 4. Section 5 
summarizes the outcome and future prospects. 

2. Stellar mass-loss and long-term galactic 
evolution 

2.1. Mass-loss rate of a coeval stellar population 

Stars of all masses lose matter via stellar winds, high- 
mass stars, in addition, in supernova explosions. Despite 
the fact that stellar winds accompany all stages of stellar 
evolution (see de Jager et al. 1988, for observations of wind 
rates across the Hertzsprung- Russell diagram), they vary 
in intensity and duration: gas is not released uniformly 
over a stellar lifetime but typically in a few relatively short 
episodes. 

The time dependence of mass-loss rate of a coeval stel- 
lar population (hereafter CSP), M{t), could in principle 
be obtained if we knew the stellar initial mass function 
(IMF), ip{mi), and the time dependence of the mass-loss 
rate of individual stars, m: 



Normalized mcss-.oss rate of o CSP 



Mass and cumulotive moss-loss o" o CS.^ 



M{t) 



m{t, rrii) ipirrii) dmi 



Since neither the former nor the latter are known with 
much precision, we shall limit ourselves to construct a 
model curve based on several simplifying assumptions ad- 
equate to the purpose of the work which is to simulate 
the effects of stellar mass- loss on the long-term evolution 
of galaxies; we do not pretend to model changes occur- 
ring on a time-scale of less than ^ 10 Myr that is short 
compared to the dynamical time at most galactocentric 
distances. 

Mass-loss for stars of different masses end evolution- 
ary stages is reviewed in more detail in Jungwiert (1998) 
where references can be found for both observations and 
models. Here, only a brief summary of the most important 
mass-loss mechanisms is given: a) in the case of low-mass 
stars (LMS, hereafter defined as stars with initial masses 
TOi < 2 M0), winds on the red giant branch (RGB) and 
asymptotic giant branch (AGB) are the most efficient; b) 
for intermediate-mass stars (IMS, 2 Mq < toi < 8 Mq), 
the AGB winds dominate; c) in the case of high-mass stars 
(HMS, TOi > 8 Mq), supernova explosions, Wolf-Rayet 
winds and main-sequence winds, contribute to the overall 
mass- loss in proportions strongly depending on m{: mass- 
loss is dominated by supernovae for ttij < 25 Mq, while 
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Fig. 1. Left: Normalized mass-loss rate, M{t)/Mi, of a co- 
eval stellar population (CSP) with Z = 0.02 and the Scalo's 
(1998) IMF. The full curve is derived from the Padua stellar 
models, the dashed curve is the fit corresponding to Eq. (2) 
with parameters from Table 1 (for r = 20 Gyr). Note that 
both scales are logarithmic. Right: Normalized stellar mass, 
Mstars/Mi (upper full curve), and the cumulative mass-loss (i.e. 
the released gas mass), Mgas/Mi (lower full curve), for the same 
CSP. The dashed curve corresponds to the normalized stellar 
mass computed for the mass-loss rate fit from the left pannel. 



for more massive stars the winds 



more efficient than 



supernovae. 

With regard to the above, we shall assume, for simplic- 
ity, that: a) for LMS all the mass- loss takes place in two 
distinct delta function events corresponding to the tips of 
RGB and AGB (since the respective winds strongly peak 
near these tips); b) for IMS and HMS, we treat the mass- 
loss as one delta function at the end of their lives: this 
reflects the AGB wind peak for intermediate-mass stars 
and SN explosion -I- winds for high-mass stars (as the 
most massive HMS, for which the winds dominate over 
supernovae, have very short lifetimes - less than 5 Myr 
for TOi > 40 M0 - the delta function representation is 
adequate) . 

Necessary ingredients to construct an approximate 
mass-loss curve of a CSP thus reduce to the knowledge of 
IMF, stellar initial- final mass relation and stellar lifetime- 
mass relation (in the case of LMS, two lifetimes are nec- 
essary, for the ends of the RGB and AGB phases) . As for 
the IMF, we use the power- law with three slopes (0.2, 1.7 
and 1.3 for initial mass ranges of 0.1-1 Mq, 1-10 Mq and 
10-100 Mq, respectively) inferred by Scalo (1998) as the 
"average" of recently published IMFs. For the initial-final 
mass relation and lifetime-mass relation, we rely on the 
Padua stellar evolutionary models (Bressan et al. 1993; 
Marigo et al. 1996) computed for the initial metallicity 
Z = 0.02. 

The mass- loss rate function, Afn(i), constructed in this 
way (the subscript "n" denotes the normalization to the 
initial mass. Mi, of a CSP) is shown in Fig. 1 (left). For 
the purpose of its efficient implementation into the N-body 
code, we carry out a least-squares fit. The power-law fit, 
Mn(t) oc t^y suggests a very simple form since it gives 
y = 1.04, very close to unity. Rather than using 1.04, 
we fix 2/ = 1, and perform a new least-squares fit with 
a hyperbola having a shift Tg from the moment tbirth, at 
which a CSP was born, to avoid the unphysical singularity 
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Table 1. Fitted mass- loss rate parameters 

r [Gyr] c,, To [Myr] 



5 
10 
15 
20 



5.55 ■ 10" 
5.47 • 10" 
5.41 • 10" 
5.35 • 10" 



5.04 
4.97 
4.92 
4.86 



Table 2. Cumulative mass- loss of a CSP 

Cumulative mass-loss at several chosen times (in Gyr) for the 
Scalo's IMF (1998) and Padua Z = 0.02 stellar models. 



Rt^o.i Rt=i Rt=b Rt=io Rt=ib Rt=2o 
0.170 0.280 0.382 0.416 0.434 0.445 



of the t ^ function at tbirth^ 



Mn{t) 



M{t) 
Mi 



Co 



t — ibirth + To 



(2) 



We restrict the least-squares fits to combinations (co, Tq) 
giving the correct total gas return (cumulative mass-loss) 
Rr over time r (see Sect. 2.2. and Table 2), corresponding 
to the time interval between the birth of a given CSP and 
the envisaged end of an N-body simulation. Results of the 
fit for T = 5, 10, 15 and 20 Gyr are summarized in Table 
1. 

The fit for the simulation length of 20 Gyr is shown 
in Fig. 1 (left) by a dashed line. It is close at all times 
and over four orders of magnitude on the M^{t) axis to 
the original dependence. It is also obvious, from Table 
1, that the fits for different r are very similar. We will 
take advantage of this when implementing the mass-loss 
into the N-body code: only one fit will be used for all the 
particles (see Sects. 3.4. and 3.5.). 

2.2. Cumulative mass-loss of a coeval stellar 
population 

While the time-dependence of mass-loss rate of a CSP is 
a vital input for our models, the cumulative mass-loss, i.e. 
the mass fraction lost by stars until time t after the birth 
of the CSP, 



Rt 



Mn{t')dt', 



(3) 



gives a useful estimate of the importance of gas recycling 
on different time-scales. 

Fig. 1 (right) shows, in the lower half of the plot, Rt of 
a CSP computed for the Scalo's (1998) IMF and the Padua 
Z = 0.02 stellar models. The curve in the upper half of 
the plot gives the fractional mass remaining in stars, i.e. 
1 - Rt. Table 2 gives Rt for t = 0.1, 1, 5, 10, 15 and 20 
Gyr. 



A discussion of mass-loss for other IMFs and stellar 
metallicities is beyond the scope of this paper. Rt for sev- 
eral other IMFs widely used in the litterature (Salpeter 
1955; Miller & Scalo 1979; Kennicutt 1983) as weU as 
for the Padua stellar models corresponding to Z=0.0004, 
0.004, 0.008, 0.02, 0.05 is tabulated in Jungwiert (1998). 
Here we note only a few facts. For Z — 0.02, the highest 
gas return over the Hubble time (let's say 15 Gyr) is ob- 
tained for the Scalo's (1998) and Kennicutt's (1983) IMFs 
(^15 Gyr — 0.43 and 0.44, respectively), followed by the 
Miller & Scalo's (1979) IMF (0.41), while the Salpeter's 
(1955) IMF gives a markedly lower value (0.30). In all the 
cases, the returned mass fraction is considerable. As for 
the relative contributions of LMS, IMS and HMS to the 
total mass-loss, they are comparable, so that no category 
can be considered as significantly dominant. For example, 
the Scalo's (1998) IMF and Z = 0.02 stellar models lead 
to ratios LMS : IMS : HMS - 1 : 1.25 : 1 for the contribu- 
tions to i?i5Gyr- The Stellar metallicity does not change 
the mass-loss of the whole CSP too much despite the fact 
that individual stars have their lifetimes and final masses 
infiucnced by Z quite noticeably (see Figs. 16 and 17b in 
Jungwiert 1998). 

2.3. Galactic gas return rate 

We define, for the purpose of simulations presented in 
Sect. 4, gas return rate (GRR) as the analog of star for- 
mation rate (SFR): GRR (units of MQyr"^) is the rate 
at which stellar mass is converted into gaseous mass. In 
general, GRR depends on the whole history of star forma- 
tion: 



GRR{t)= SFR{t') Mn{t~t')At', 



(4) 



where M^ {t — t') is the (normalized) mass-loss rate, at 
time t, of a stellar population born at t' . Note that when 
GRR refers to a specified area, SFR in Eq. (4) does not 
refer to the same area: GRR at time t is not determined 
by the star formation history of this area but by the star 
formation history of stellar matter that is now present in 
this area but has, in general, various orbital histories. 

3. N-body code with star formation and 
time-dependent stellar mass-loss 

3.1. Gravitational interaction and equations of motion 

We use the standard particle-mesh (PM) technique 
(Hockney 1970; Hohl 1971; for a review, see: Hockney 
& Eastwood 1981; Selwood 1987) to compute the grav- 
itational potential on a Cartesian grid via Fast Fourier 
Transform (FFT). Forces at grid points are then evalu- 
ated by differencing the potential. Both the assignment of 
particle masses to grid points and the finding of forces be- 
tween them are carried out by the bilinear "cloud-in-cell 
(CIC)" interpolation. The equations of motion arc inte- 
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grated by means of the leap-frog algorithm (Hockney & 
Eastwood, 1981). 

The gravitational softening that reduces interparticle 
forces at short distances (see Pfenniger & Friedli 1993) 
is that of a Plummer sphere. The choice of the softening 
length is done with regard to the criteria laid down by 
Romeo (1994, 1997) to optimize the stability and relax- 
ation properties of N-body disks. 

3.2. Particle species 

There are four particle species in our code. Apart from 
relatively massive stellar particles (hereafter referred to 
as "standard stars"), typically representing millions of 
stars, and relatively massive gaseous particles (hereafter 
"clouds") with characteristic masses of 10'"' — 10® Mq, we 
introduce low-mass stellar and low-mass gaseous particles 
(hereafter nicknamed "starlets" and "cloudlets" ) that me- 
diate the exchange of matter - converted from gaseous to 
stellar phase and vice versa by star formation and stel- 
lar mass-loss - between the standard stars and clouds. 
The definition of mass boundaries between cloudlets and 
ordinary clouds and between starlets and ordinary stars, 
and thereby the mass ranges for cloudlets and starlets, are 
precised in Section 4.1. (item 13). 

The standard clouds are meant to simulate the cold 
gas component. They are modelled as finite-size particles 
that undergo physical collisions in which they dissipate 
a fraction of kinetic energy. The relative velocity of two 
colliding clouds is reduced by factors /3r, A; the former 
multiplies the velocity component along the line joining 
the two clouds (radial component) , the latter perpendicu- 
lar to it (tangential component) (for pioneering use of this 
"sticky-particles" technique in simulations of galaxies, see 
Schwarz 1981; among more recent applications and dis- 
cussion on values of /3i., /3t, see e.g. Jungwiert & Palous 
1996). 

Star formation is allowed to proceed in our cold gas 
component in the way described in Sect. 3.3. It results in 
the production of the above mentioned starlets. These in 
turn join the standard stellar component by merging with 
a standard star provided some is near by (typically 200 
pc; to be precised below). On the other hand, standard 
stars lose continuously mass in the form of the cloudlets - 
as detailed in Sect. 3.4. - that join the cold gas component 
by merging with a standard cloud provided some is near 
by. The condition of particle "nearness" for star-starlet 
or cloud-cloudlet merging ensures that the mass exchange 
between the stellar and gaseous components is local. 

In usual situations the existence of our starlets and 
cloudlets will be ephemeral since they will merge with a 
nearby standard particle immediately after they are born. 
This is technically equivalent to a direct exchange of a 
fraction of mass between a standard stellar and a stan- 
dard gaseous particle. However, there are important and 
frequent astrophysical situations in which stars and cold 
gas are segregated spatially. For example, in some areas 



cold gas can be depleted by intense star formation or swept 
away by gravity torques. There can also be regions of cold 
gas without any or with only small amounts of stars, for 
instance the whole galaxies in early evolutionary stages 
or, in the case of present-day galaxies, the area beyond 
the periphery of optical disks. 

Our starlets/cloudlets enable to follow such situations. 
If a starlet/cloudlet is produced in an area (defined be- 
low) with no standard star/cloud, it leads an independent 
and possibly long-term existence until it meets the respec- 
tive standard particle (meanwhile, the starlet is subject 
to mass- loss as standard stars and it produces cloudlets). 
The long-living starlets/cloudlets can also meet another 
starlet/cloudlet. Such events are treated by starlet-starlet 
and cloudlet-cloudlet mergings. In the course of time, 
a starlet/cloudlet grown by successive mergings of star- 
lets/cloudlets might approach the mass typical of stan- 
dard stars/clouds. If this occurs, the starlet/cloudlet is 
converted into a standard star/cloud. Conversely, a stan- 
dard cloud is changed over to a cloudlet if its mass falls, 
due to the starlet production, below the value typical of 
standard clouds. In contrast, we do not consider a con- 
version of standard stars into starlets since the mass of 
a standard star remains of the same order during a sim- 
ulation: it cannot decrease below the quantity 1 — i?t^„j 
that is larger than ~ 0.55 for any simulation length tend 
shorter than 20 Gyr (see Table 2). 

In this respect, the cloudlets can be viewed as a very 
crude representation of the warm/hot galactic gas: they 
are products of star formation and do not participate in 
the star forming process unless they "cool" in areas of high 
gas density (either by merging with a standard cloud or 
with many other cloudlets). 

In summary, the four particle species have the follow- 
ing properties. Standard stars produce cloudlets and can 
absorb starlets. Starlets can be absorbed by standard stars 
or merge with other starlets; they produce cloudlets; if 
their mass exceeds a given value (Mgs, ~ 10^ — 10® Mq, 
see Sect. 4.1., item 13), they are converted into standard 
stars. Standard clouds undergo mutual, partly inelastic 
collisions; they produce starlets and can absorb cloudlets; 
if their mass falls below a given value (Mcc, ~ 10^ Mq, 
Sect. 4.1., item 13), they are converted into cloudlets. 
Cloudlets can be absorbed by standard clouds or merge 
with other cloudlets; if their mass exceeds the above men- 
tioned Mcc, they are converted into standard clouds. 

All the particle productions, absorptions, mergings and 
collisions happen in the mass and momentum conserving 
manner. 

Perhaps the most questionable interaction in our code 
is the merging of a starlet with a standard star or another 
starlet. The real stellar fluid cannot dissipate energy; nev- 
ertheless, if new stars are born with a low velocity dis- 
persion (as expected from the kinematics of parent gas 
clouds), the velocity dispersion of the real stellar fluid as 
a whole can decrease. We merge the particles for computa- 
tional reasons, so as to keep their number low (see below); 
we are unable to follow every new born stellar particle 
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independently. In this context, the merging of standard 
stars with starlets can be viewed as very crudely mim- 
icking the rejuvenation of the stellar fluid by newly born 
stars and the related decrese of its velocity dispersion. 
The merging of one starlet with another can seem more 
problematic but this process is relatively rare in simula- 
tions we present. However, we admit that the merging of 
stellar particles could in general bias the evolution of the 
stellar orbital structure, especially if the merged particles 
had very different velocities. This problem certainly merits 
further investigation. 

As indicated above, the interparticle interactions (col- 
lisions, mergings) require the ability to identify nearby 
particles. The same is true for the star formation process 
described in the subsequent subsection. To find nearby 
particles in a CPU-efficient way, we divide a galaxy into 
a 2D rectangular grid (hereafter "interaction grid" ; for 
its specifications see Sect. 4.1., items 11 and 6); four lists 
(standard stars list, starlets list, standard clouds list and 
cloudlets list) are associated with it: within each of them, 
the particles are continuously sorted in such a way that 
their indices are stored in the order of the increasing index 
of the grid boxes. This permits a fast identification of all 
the particles in a given box at every time step. Particles 
located in a given box are considered as "nearby" in the 
above used sense. It is true that two particles from distinct 
neighbouring boxes may be closer to each other than par- 
ticles within a same box. However, our concern is basically 
to have at most one starlet and one cloudlet per box and 
to achieve this in a computationally simple manner. The 
fact of ignoring interactions behind the box boundary is 
dynamically unimportant since it happens on distances 
below the dynamical resolution of the code. 

The star-starlet, starlet-starlet, cloud-cloudlet and 
cloudlet-cloudlet mergings permit to keep a low number 
of particles. If starlets/cloudlets were not merged, they 
would soon proliferate in an inacceptable number since 
at each time step at which the star formation and mass- 
loss procedures are applied, roughly N^oxes {^ 10^ — 10^) 
starlets/cloudlets are created. The merging rule usually 
leads to mergings of nearly all starlets/cloudlets with stan- 
dard stars/clouds within the time step in which these star- 
lets/cloudlets were created (see Sect. 4.3.) except areas or 
evolutionary stages where or when there is not enough 
stars or cold gas. 

3.3. Star formation scheme 

We impose a generalization of the Schmidt law for star for- 
mation (Schmidt 1959; for a recent review, see Kennicutt 
1998), i.e. the surface density of star formation rate, fJ^sFR, 
is proportional to a power of the gas surface density: 

fJ'SFR = ^A^gas- 

Kennicutt (1998) gives a = 1.4 as the "best-fit" for 
the obscrvationally determined Schmidt law relating the 
mean surface density of SFR to the mean surface density 
of cold (molecular + atomic) gas in both galactic disks 



and circumnuclcar regions of disk galaxies. In our model, 
we keep the freedom of choosing a as an input parameter 
(typically 1 or 2). The Schmidt law is applied with a time 
step AtsF (longer than the time step used to integrate the 
equations of motion; Sect. 4.1., item 12) separately in each 
of the boxes of the interaction grid. 

Inside a given box, star formation proceeds as follows: 
1) The standard clouds list is used to compute the cold 
gas surface density. 2) If this density is non-zero (i.e. there 
is at least one standard cloud), a starlet is created with 
the mass computed according to the Schmidt law. 3) The 
starlet mass is subtracted from the standard clouds in 
amounts proportional to their masses. 4) The position and 
velocity of the starlet are found in such a way that the to- 
tal momentum of the standard clouds and of the starlet 
is conserved. 5) Standard stars, if any, are found using 
their list. If there are some, one of them is randomly cho- 
sen to merge with the starlet. If there is no standard star, 
the starlets list is searched through for present starlets, 
if any, that immigrated from neighbouring boxes. If the 
search is positive, all the starlets in the box merge into 
one. Otherwise the new starlet will temporarily lead an 
independent existence as described in the previous sub- 
section. 

3.4. Stellar mass- loss scheme 

In Sect. 2 we have evaluated the normalized mass-loss rate, 
Mn(t), of a coeval stellar population. Here our goal is to 
implement the mass-loss into the N-body code. 

As described in two previous subsections, new stellar 
mass forms, in our code, from gas as starlets that can 
merge with nearby standard stars or other starlets. Our 
typical stellar particle (either a standard one or a starlet) 
therefore docs not represent a CSP but has a complex star 
formation history since it consists of mass of starlets born 
at various times in addition to the mass it was assigned 
in the beginning of the simulation. 

The exact treatment of the mass-loss - for the code 
purposes discretized in time with a step A^ml (Sect. 4.1., 
item 12) - would then require to keep in memory, for each 
stellar particle, its star formation history, i.e. initial masses 
of starlets that were incorporated into it and the times at 
which these starlets were born. The mass lost by stellar 
particle J at time t; due to incorporated starlets would 
read: 



AM,/ (to = 



'ml k<l 



m,j{tk)-Mu{t' ~tk)dt' = 



= ^m,j{tk) ■ AMniti 



-tk) 



(5) 



where m<,j{tk) is the mass with which a starlet, incor- 
porated into J aX t < ti, was created at time tk (the in- 
corporation could have happened later than at tk), and 
AMn(ti — tk) is the fraction of mass lost over ti — Atuh 
by the stellar mass born at tk- Computing this sum for 
every stellar particle (iVstars+stariots ^ 10'"') in each step 
(-^stops ^ 1000 if mass-loss is applied every A^ml ^ 10 
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Myr over the Hubble time) and keeping in memory the 
arrays of to* and tk would be computationally very ex- 
pensive. 

We therefore introduce a much simpler manner to eval- 
uate the mass-loss of stellar particles getting benefit from 
the fact that Mn(i) was shown to be reasonably well fitted 
by a simple analytical function with only two parameters 
(Eq. 2). The trick consist in using the same functional form 
for all stellar particles disregarding their star formation 
histories, while individually adjusting the two parameters 
whenever a new starlet is incorporated. 

It is thus necessary to define how to change the rel- 
evant parameters when two stellar particles of different 
ages and masses merge, and to check whether or not such 
a change plausibly represents the sum of their mass-loss 
rates. Formulated mathematically, we want to find, for a 
particle J formed at time tj by merging of particles R and 
S, born at t^ and ig, parameters Cj and Tj such that the 
mass-loss rate of J (hereafter called the combined mass- 
loss rate), expressed as 



Mj{t) 



Cj 



t-tj+ Tj ' 



(6) 



would have similar time dependence and total gas return 
as if particles R and S were followed independently. Eq. 
(5) is then replaced by the corresponding integral of eq. 
(6): 

AMj(t,) = ~Cj ln(l - -^^^) (r) 

ti — tj + 1 J 

Since there are two free parameters, Cj and Tj, we can 
impose two conditions on the combined mass-loss. After 
trying several possibilities (for a more detailed analysis, 
see Jungwiert (1998), pages 52-58) we require Eq. (7) to 
give the correct cumulative mass-loss at two time mo- 
ments: at the end of the simulation, tend, and at an in- 
termediate time, tim, between tj and tend; the choice of 
tim = tj + 0.2 • (tend — tj) proved to give good results. 
These two conditions can be expressed as 



Mj(t)df : 



tj 



MR{t)dt+ / Ms(t)df, (8) 



where tx denotes either tend or tim. 

A priori, it was not ensured that our combined mass- 
loss approximation well represented the real sum of par- 
ticles' mass-losses. However, the tests presented in Sect. 
3.6. verify the adequacy of our approach. 

Knowing how to combine the mass-loss rates of indi- 
vidual particles when they merge, it remains to be said 
what values of Cj and Tj we assign to particles that do 
not form by merging, i.e. to starlets at the moment of their 
birth and to standard stars, if any, present in the begin- 
ning of the simulation. Such can be viewed as representing 
coeval stellar populations. Comparing eqs. (2) and (6), we 
see that Tj is simply Tq and Cj is cq multiplied by the 
mass of a particle at the moment of its birth, tbirth = tj. 
Values of cq and Tq are chosen from Table 1 (as we have 
mentioned in Sect. 2.1., the fits presented in Table 1 are 
very close to each other; we will use that for t = 10 Gyr). 



For the standard stars present in the beginning of the 
simulation, there is an additional issue consisting in defin- 
ing their tbirth! we postpone it to a special subsection be- 
low. 

The last point is to explain what happens with masses 
AM J lost by individual stellar particles at every time step 
AtML- We use a strategy similar to that for star formation 
(Sect. 3.3.). The procedure runs through the interaction 
grid box by box: 1) All the standard stars and starlets are 
found using the respective lists and their AMj (computed 
according to Eq. (7)) are summed up. 2) If the sum is 
non-zero (i.e. there is at least one stellar particle in the 
box), the mass is attributed to a new cloudlet. 3) It is 
positioned to the center of mass of contributions AMj 
coming from individual stellar particles and its velocity 
is computed in such a way that the total momentum of 
stellar particles in the box and of the cloudlet itself is 
conserved. 4) The standard clouds, if any, are found using 
the clouds list. If there are some, one of them is randomly 
chosen to merge with the cloudlet. If there is no standard 
cloud, the cloudlets list is searched through for cloudlets, 
if any, that immigrated from neighbouring boxes. If the 
search is positive, all the cloudlets in the box merge into 
one. Otherwise the new cloudlet will temporarily lead an 
independent existence. 



3.5. Mass-loss of pre-existing stars 

The purpose of this subsection is two-fold: first, to drive 
attention to the fact that when a simulation of a galaxy 
starts from a state in which a stellar component is already 
present, as is the case of numerous published works, it is 
natural and important to take into account its mass-loss; 
second, to define tbirth and Cj for standard stellar parti- 
cles representing this component (hereafter "pre-existing 
stars" and "pre-existing stellar component"). 

To know at what rate the pre-existing stellar compo- 
nent should lose mass during the simulation, an assump- 
tion must be made about the star formation history of 
the model galaxy. Here we restrict our attention to two 
extreme cases, a constant star formation rate over some 
epoch prior to the start of the simulation (defined by 
t = 0) and a starburst at some moment prior to the start 
of the simulation. To be more specific, we suppose a sim- 
ulation of length tend = 10 Gyr and calculate what gas 
mass, AMgas, is released during the simulation (i.e. be- 
tween t = and tend) by the pre-existing stellar compo- 
nent. Table 3 gives this mass (relative to the mass of the 
pre-existing stellar population at t = 0) for two constant 
SFR histories (lasting 5 and 10 years) and two starbursts 
(occuring 5 and 10 years prior to the start of the simula- 
tion): 

The values of AMgas range between 5 and 15% of 
the pre-existing stellar mass. This can represent a huge 
amount of gaseous matter in absolute terms. For exam- 
ple, starting a simulation with a stellar component of 
4 • 10^" Mq (as will be the case of our illustrative sim- 



B. Jungwiert et al.: Continuous stellar mass-loss in N-body models of galaxies 



AMgaB(tc„d = 10Gyr) 


E ° ■ 

* o r 


Matars (t=0) 

0.153 
0.115 


0.084 


Y '■ 


0.050 


5 1 



Table 3. Gas released by pre-existing stars 

pre-simulation SF; 

const. SFR (-5 to Gyr) 
const. SFR (-10 to Gyr) 
starburst (-5 Gyr) 
starburst (-10 Gyr) 



ulation in Sect. 4) implies that the gas mass released from 
these stars between t — and tend = 10 Gyr would range 
between 2 • 10^ Mq and 6 • 10^ Mq which is typical of gas 
mass in disks of observed present-day spiral galaxies. The 
mass-loss of pre-existing stars is thus quite important and 
should not be neglected in simulations. 

Technically, to be able to use Eq. (7), we need to assign 
to each of the pre-existing star (having mass Mq at i = 0) a 
birth time, ibirth = tj, and a birth mass, Mbirth (the latter 
will give us Cj after multiplying by cq). In the starburst 
case, all such particles have ibirth simply equal to the time 
moment of the starburst; in the case of a constant SFR 
history, we draw tbirth for each pre-existing particle ran- 
domly within the period of constant star formation epoch. 
Once we assign ibirth? -^birth is found by integrating Eq. 
(2) backwards in time (from i = to ibirth < 0). 

Pre-existing stars are thus supposed to have already 
lost mass Mbirth — Mq prior to the start of the simula- 
tion, however their mass-loss will continue also during the 
simulation. 



3.6. Tests 

Several tests of our mass-loss algorithm are described in 
Jungwiert (1998). Especially, it has been tested that the 
combined mass-loss (Sect. 3.4.) of a stellar particle built 
by merging ("composite particle" for the purpose of this 
subsection) of two or more stellar particles ( "merged parti- 
cles" ) well approximates the sum of individual mass- losses 
for many mass and age ratios of merged particles and for 
many subsequent mergings. Five examples are shown in 
Figs. 2-4. In each row, the left pannel shows the mass-loss 
rate of a composite particle while the right pannel shows 
its mass and cumulative mass-loss (i.e. the gas mass it has 
released). Full lines represent quantities computed as if 
all the merged particles were followed independently (i.e. 
according to Eq. 5) while the dashed hues correspond to 
the combined mass-loss approximation (eqs. 6-8) used in 
our code. All the quantities are normalized to the sum of 
initial masses of merged particles. 

Fig. 2 corresponds to a merging of two stellar particles 
of the same initial mass. The first one is born at t = 0, 
the second at t — 5 Gyr when the merging takes place. 
Fig. 3 shows the case of multiple mergings, equidistant 
in time, with all the merged particles having the same 
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Fig. 2. Merging of two stellar particles born at i = and 
5 Gyr. Left: mass-loss rate of the composite particle; right: 
mass (Mstars, upper curve) and cumulative mass-loss (Mgas, 
lower curve) of the composite particle. Full lines in Figs. 2-4 
correspond to exact sums of mass-losses of individual particles, 
dashed lines represent our combined mass-loss approximation. 




Fig. 3. Like Fig. 2 but for multiple mergings of particles hav- 
ing the same initial mass and born equidistantly in time. The 



particles merge at their birth. Top: A'n: 



= 10; bottom: 



A^morgc = 10000. 



initial mass. The number of particles, TV, 



merge 5 



merged into the composite particle is 10 in the top row 



and 10 GOO in the bottom row. Fig. 4 also corresponds to 
A^merge = 10 (top row) and 10 000 (bottom row) but now 
the first particle is ten times, respectively ten thousand 
times more massive than the others it merges with. In all 
the examples presented in Figs. 2-4 the combined mass- 
loss approximation well represents the exact mass-loss. A 
few additional examples can be found in Jungwiert (1998, 
pages 52-58). 

We have also carried out several runs of the N-body 
model in which the motion of particles was inhibited 
while the star formation and mass-loss schemes, producing 
cloudlets and starlets and combining them with standard 
clouds and stars, were operational. We have verified that 
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Fig. 4. Like Fig. 3 (i.e multiple mergings of particles born 
equidistantly in time) but the initial mass of the first particle, 
mi,i, is higher than initial masses, mi^othors, of other particles it 
merges with. Top: A'mcrgo = 10, rrnA — 10 m-i, others; bottom: 



Nrr 



= 10000, mi.i = 10000 mi.othors. 



these runs yield temporal evolutions for SFR and GRR 
close to numerical solutions of models not using any par- 
ticles. 



4. An example: a disk galaxy 

As an illustrative application of our code, we present 4 
simple 2-D models of a disk galaxy having no gas in the 
beginning of simulations. The first two models (denoted 
A and B) are rather artificial but useful reference cases: in 
model A stellar mass-loss is supressed, the disk remains 
100% stellar (classical purely stellar N-body model): in 
model B mass-loss from initial stellar particles produces 
gas but star formation is supressed. The standard mod- 
els (denoted CI and C2) are those in which both star 
formation and stellar mass-loss operate: model CI uses 
the linear Schmidt law {a ~ 1), model C2 the quadratic 
one (a = 2); stellar mass- loss is due both to initial and 
newly formed stellar particles. Table 4 summarizes basic 
attributes of the models. The simulation length is 10 Gyr 
for all of them. 



4.1. Input data 

For the sake of completness and reproducibility, we enu- 
merate, before describing the simulations, the code input 
parameters and specify their numerical values. The ma- 
jority of them will be common to a larger set of models 
that will be published in a separate paper. Items 1-7 refer 
to parameters characteristic of collisionless N-body mod- 



Basic attributes of models 






star formation 


mass-loss 


IVIodel 






A 


no 


no 


B 


no 


yes 


CI 


yes, a — 1 


yes 


C2 


yes, a — 2 


yes 



els, items 8-13 are specific to our gas dynamics and star 
formation/mass-loss schemes. 

1) Initial stellar surface density (2D Toomre-Kuzmin 
disk: Mdisk = 4-10^° Mq, scale-length of 4 kpc, truncation 
radius of 16 kpc); 

2) Initial number of stellar particles, A'stars — 50000; all 
of them have initially the same mass Mo — Mdisk/A^stars = 
8-105 IVIq. 

3) Initial velocity distribution. We fix the Toomre 
(1964) stability parameter Qstars — 1-5 (constant in 
space), and then obtain velocities of particles by finding 
the time-independent axisymmetric solution of the Jeans 
equations; the ratio of azimuthal-to-radial velocity disper- 
sions is computed using the epicyclic approximation (see 
e.g. Binney & Tremaine 1987); 

4) Bulge and halo masses (Mbuigo = Afdisk, Afhaio ~ 5 • 
10^^ M0) and scale-lengths (Plummer spheres of 1 and 25 
kpc, respectively). The bulge is implemented essentially to 
partly stabilize the inner disk, the halo to keep the rotation 
curve fiat in the outer disk. Both components (introduced 
by analytical formulae) are treated as unevolving; 

5) Gravitational softening length (Sect. 3.1.), s k, 0.2 
kpc. It corresponds to the Romeo's (1994, 1997) "safe" 
softening computed for the center of our disk. 

6) Number of divisions of the Cartesian grid used to 
compute the gravitational potential (Sect. 3.1.): ND = 
256 (in each dimension). This number relates to the "ac- 
tive" grid (the computation proceeds on a grid doubled in 
each dimension as required by the FFT technique). The 
linear size of a grid cell, Agrid , is set equal to the softening 
length. The above choice of ND and s gives the linear size 
of the active grid, -Dgrid = 50 kpc (allowing for maximum 
disk radius of 25 kpc) ; 

7) Time step for the integration of equations of motion, 
Af . It is constant over the whole disk. Its value varies with 
time in such a way that it is 1/25 of the minimum orbital 
period. Initially At = 0.6 Myr. It decreases as the central 
mass concentration increases; 

8) Parameters for cloud collisions (Sect. 3.2.): a) /3t = 
1.0. This choice assures the conservation of angular mo- 
mentum in a collision; b) (3-^ = 0.6. This choice is rather 
arbitrary. Varying /Sj. falls beyond the scope of this pa- 
per. We just note that a positive sign for /J,, means, in our 
notation, that the direction of the cloud-to-cloud momen- 
tum in a collision is not reversed, unlike in most works 
implementing this coUisional scheme. For more comments 
on values of /J,., see Jugwiert & Palous (1996); 
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Fig. 5. Disk evolution in the purely stellar model A. In Figs. 5-7: numbers in the top right corner give time in Myr; only 20% 
of standard stellar particles, i.e. 10000, arc shown for clearness. 
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Fig. 6. Disk evolution in the standard model CI. Top: standard stars; middle: standard clouds; bottom: standard clouds 
zoom to the central region. 



9) Two parameters for the Schmidt star formation law 
(Sect. 3.3.): we choose a = 1 or 2 and compute the pro- 
portionality constant A in such a way that the SFR value 
coincides with that of the Kennicutt's (1998) fit (a = 1.4) 
at cold gas surface density of 5 M© pc~^ (a value typical 
of present day-spiral galaxies) ; 

10) Parameters related to the stellar mass-loss rate 
(Sects. 3.4. and 3.5): to assign ibirth to initial stellar parti- 
cles, we use a simplifying arbitrary assumption that SFR 



was constant over a period Tgp = 5 Gyr before the be- 
ginning of the simulation, cq and Tq (that determine Cj 
and Tj of individual standard stars and starlets) are taken 
from Table 1 (for t = tend = 10 Gyr); 



11) The interaction grid (Sect. 3.2.) used for cloud col- 
lision, star formation and stellar mass-loss schemes has 
the same size and the same number of cells as the grid 
described in item 6; 
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12) Time steps AtsF (Sect. 3.3.) and AiML (Sect. 3.4.) 
for applying star formation and mass-loss schemes. We fix 
AtsF — A^ML = 6 Myr (10 times the initial integration 
step At); 

13) Masses M^s and Mcc for star-starlet and cloud- 
cloudlet conversions (Sect. 3.2.). We fix Mss = (1 — 
-Rt=ioGyr) • Mo ~ 3.3 • 10^ M©, i.e. equal to the mini- 
mum mass to which the cumulative mass-loss Rt (Table 
2) can drive a stellar particle with initial mass Mg. As for 
Mcc, we set 10^ Mq; the overwhelming part of cold gas in 
galaxies seems to dwell in clouds with higher masses than 
this value. 



4.2. Results 

The disk evolution in the purely stellar model A and in the 
standard model CI is shown in Figs. 5 and 6, respectively. 
Stellar disks (Fig. 5 and upper row of Fig. 6), initially ax- 
isymmetric, first develop transient spiral arms and then 
a long-lasting bar. Meanwhile, in model CI, stars release 
gas in the form of cloudlets. These gradually merge and in 
regions of high gas density form standard clouds (see Sect. 
3.2.) that kinematically react to the non-axisymmetric 
potential and accumulate predominantly in a dense nu- 
clear ring (Fig. 6, bottom row) and in an outer ring (Fig. 
6, middle row), connected respectively to the inner and 
outer Lindblad resonances (ILR, OLR). A hint of another, 
weaker gaseous ring, surrounding the stellar bar and re- 
lated to the 4/1 resonance, can also be seen. 

For models B and C2 (Fig. 7, left and right pannels, re- 
spectively) , we show only an advanced evolutionary stage 
(i = 9 Gyr) that can be compared with the last column 
snapshots of Figs. 5 and 6. In model B, one can note that 
the bar is shorter and rounder (than in all the other mod- 
els), the outer ring is less regular (than in CI and C2 
models) while the nuclear ring has collapsed to a smaller 
size. All these differences are related to a huge gas release 
and associated bar weakening (for a review, see Combes 
2001) that are discussed below. The disk structure and 
ring shapes in model C2 are rather similar to those of 
model CI, except that the standard (i.e. relatively mas- 
sive) clouds in the nuclear ring are considerably depleted 
by a more intense star formation (resulting from a higher 
star formation efficiency at high gas densities). 

As indicated under item 10 of Sect. 4.1., we have cho- 
sen the assumption of constant SFR over 5 Gyr preceding 
the start of the simulation. Given the initial disk stellar 
mass of 4 • 10^" M©, Table 3 implies that the gas mass 
released by the pre-existing stars over the simulation is 
about 6 • 10^ Mq. While our choice of the star formation 
history is arbitrary and favouring a rather huge mass-loss. 
Table 3 indicates that other choices would lead to a gas 
restitution of the same order. 

We show in Fig. 8, separately for the central disk re- 
gion (here defined by i? < 2 kpc) and for the outer disk 
(i? > 2 kpc), temporal evolutions of SFR, GRR, total 
(i.e. stellar + gaseous mass) and gas mass fraction. The 
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Fig. 7. Models B (left pannels) and C2 (right pannels) at i = 9 
Gyr. Top: standard stars; middle: standard clouds; bottom: 
standard clouds - zoom to the central region. 



gas restitution by pre-existing stars in models B, CI and 
C2 has important effects. First, it opens the way to star 
formation. Second, the presence of gas affects the radial 
transport of matter. 

As can be seen in Fig. 8, SFR is of the order of 
1 Mq yr~^ in the disk as a whole and, in the outer disk, 
it is pretty constant over the whole simulation, except 
early stages. The relative constancy of SFR is a result 
of the time-dependent recycling as was demonstrated by 
Jungwiert (1998) on a larger set of models, starting from 
either a two-component disk with Mgas/Aigtars = 0.1 or 
from a purely gaseous disk. As for GRR, we note that, 
except early stages, contributions of the pre-existing stel- 
lar population and of stars formed during the simulation 
are comparable (the latter contribution corresponds to the 
difference between GRR curves for models CI (or C2) and 
for the non-SF model B). 

The transport of matter towards the center is quanti- 
fied in Table 5 that shows the mean transfer rate of total 
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Fig. 8. Models A, B, CI, C2 - left: outer disk {R > 2 kpc), 
right: central region (R < 2 kpc); top: GRR (thin lines) and 
SFR (thick lines), middle: Mgas/Mtot, bottom: Mtot- 




tirne [Gyr] 

Fig. 9. Bar strength Qt. 
models CI and C2. 



tirne [Gyr] 

left: models A, B and CI; right: 



(stellar + gaseous) mass into the central region (i? < 2 
kpc) and the relative strength of the bar, Qtot (defined 
as the ratio of the maximum bar tangential force over 
azimuthally averaged radial force of the total force field; 
Combes & Sanders 1981; Buta & Block 2001) at 10 Gyr. 
The radial tranport is the least efficient in model A, 
as expected, since there is no gas. Model B is the most 
gas rich since mass-loss operates but star formation is 



Table 5. Mean mass transfer rate (into R < 2 kpc) and 
bar strength 



dAfto 



>inflow [ Moyr ^] Qtot (10 Gyr) 



0.30 
0.14 
0.27 
0.30 



Table 6. Gas mass fractions at 10 Gyr 

Afga3(-R<2kpc) Mga3(-R<2kpc) 



[odel 


-- dt 


^^mnow 


A 




0.17 


B 




0.27 


CI 




0.33 


C2 




0.34 



[odel 


Mtot (fl<2 kpc) 


Mga3(fl<20kpc) 


B 


0.260 


0.46 


CI 


0.055 


0.45 


C2 


0.006 


0.11 



inhibited, however it is not the one with the largest in- 
flow. The most efficient are models CI and C2. This be- 
haviour can be understood with the help of Fig. 9 that 
shows evolutions of the bar strength Qtot- In model A, 
the bar instability appears after about 1.5 Gyr, the bar 
strength steeply grows until about 3 Gyr, then saturates 
at the value of ^ 0.33 and after 5 Gyr enters a phase 
of only slow secular weakening (at 10 Gyr, Qtot = 0.30). 
The bar evolution in models B and C differs from model A 
in two respects. First, the bar appears earlier (instability 
starts before t — 1 Gyr and the saturation is reached at 
^ 2 Gyr) since the gas production helps to destabilize the 
disk. Second, after saturation, there is more important 
weakening of the bar, especially pronounced in the gas 
rich model B (bar strength decreases from 0.35 to 0.14). 
This is connected to the known fact that a too impor- 
tant accumulation of matter near the center disturbs the 
bar orbital structure (Pfenniger & Norman 1990; Hasan & 
Norman 1990; Friedli & Benz 1993). The high gas fraction 
in the disk of model B thus first leads to a rapidly growing 
central concentration of gas (in the nuclear ring) which in 
turn disturbs the bar thus slowing down the future mass 
transfer (see Heller & Shlosman 1996). 

The recycling models CI and C2 thus provide, from 
the point of view of the efficiency in transferring matter 
inwards, the best combination of gas mass fraction and 
bar force, both quantities being interdependent. 

Table 6 gives two gas mass fractions at the end of sim- 
ulations: 1) ratio of gas-to-total mass inside the central 
region (here, "total" means gaseous -I- stellar but without 
the bulge contribution); 2) ratio of gas masses inside the 
central region and in the disk as a whole. In model B, 
about one quarter of the central region dynamical mass 
is in the form of gas. Note also that in models B and CI 
roughly one half of the gas mass is located in the central 
region, i.e., more or less in the nuclear ring. 
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4.3. Computational efficiency and statistics of particle 
species 

The presented algorithm for the continuous stellar mass- 
loss has a low computational cost. In the standard model 
CI, with 50 000 standard stars and active grid of size 256 x 
256, the mass-loss and star formation related subroutines 
used about 12% and 3%, respectively, of the total CPU, its 
main part being spent on calculating the gravitational po- 
tential and advancing the particles. The simulations were 
carried out on a relatively slow one-processor PC Pentium 
II 366 MHz with one time step, Ai, taking on the average 
0.9 sec. 
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Fig. 10. Masses (left) and numbers (right) of different particle 
species (for convenience, the mass and number of standard 
stars are scaled down by factors of 20 and 3, respectively, while 
masses of starlets and cloudlets are scaled up by a factor of 10. 
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in standard stars slowly varies due to their mass-loss and 
absorption of starlets. 

The number of standard clouds quickly (on a time- 
scale of ~ 500 Myr) grows (due to the initially copious 
stellar mass-loss) from to ^ 2000 and then continues 
climbing only mildly (to 3300 at t = 10 Gyr). On the 
other hand, the overall mass in these particles reaches a 
maximum {r-^ 1.8 • 10^ M0) at about 4 Gyr and then 
declines (to ^ 1.3 • 10^ M© at i = 10 Gyr) due to intense 
star formation and decreasing mass-loss (cf. Fig. 8). 

The number of starlets and the total mass in them are 
very low throughout the simulation (^ 200 and 10® M0) 
and negligible compared to standard stars. The starlets 
arc thus almost virtual particles: despite a big number 
of them (of the order of 1000; one new starlet per box 
which hosts at least one standard cloud) is formed every 
time step A^ml (Sect. 4.1., item 12), the overwhelming 
majority is immediately incorporated to nearby standard 
stars. The only surviving starlets are located in outer parts 
of the galactic disk (Fig. 11, left) where the number of 
standard stars is low for the star-starlet mcrgings to work 
with the 100% efficiency. 

The situation is somewhat different for cloudlets since 
their number is relatively high, between 10 000 and 13 
000. However, their total mass is (except very early evo- 
lutionary stages) much lower (it is scaled up by a fac- 
tor of 10 in Fig. 10 !) than the mass in standard clouds. 
Initially, the cloudlets happen to gather ~ 10® Mq due 
to the mass- loss from initial standard stars, however they 
efficiently lose it in favour of the standard clouds (initially 
only by cloudlet-cloud conversions, later on also by cloud- 
cloudlet mergings). Cloudlets survive only in regions of 
low total gas density and most of the time harbour only a 
few 10^ M0. Looking at Fig. 11 (right panel), one sees that 
the distribution of cloudlets is largely complementary to 
that of the standard clouds: cloudlets are depleted in high- 
gas-density regions, i.e. in the nuclear and outer gaseous 
rings. Note also that the cloudlet distribution extends be- 
yond the outer edge of the standard cloud distribution (cf. 
Fig. 6). 

Overall, our star-formation and mass-loss schemes 
with four particle species keep the total number of all par- 
ticles (~ 65 000) close to the initial value (50 000) and the 
various types of the particle species interactions do not 
put a heavy load on the CPU time. 



Fig. 11. Distribution of starlets (left) and cloudlets (right) in 
model CI. 

Fig. 10 gives the statistics (masses and numbers) of 
the four particle species for the model CI. Fig. 11 shows 
spatial distributions of starlets and cloudlets. 

The number of standard stars was constant. In princi- 
ple, it cannot decrease since these particles have no way to 
disappear (except escapes that however did not occur) but 
it could rise due to starlets-star conversion after a starlet 
reaches the critical mass Mgs (Sects. 3.2. and 4.1., item 
13); such conversions however did not happen. The mass 



5. Concluding remarks 

We have presented a gas dynamical N-body model includ- 
ing star formation and time-dependent stellar mass- loss. 
The implementation of the latter feature is, both physi- 
cally and technically, the most innovative element of our 
simulations and, as such, is the major result of this work. 
The related description should allow other people study- 
ing the galactic evolution by means of N-body simulations 
to transpose our mass-loss scheme into their codes. 

Our code is obviously not free of shortcomings, for in- 
stance the ignorance of gas cooling and heating, inher- 
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ent to sticky-particles representation of gas (we however 
recall - see Sect. 3.2. - that our subdivision of gas parti- 
cles into standard clouds and cloudlets takes into account, 
on a very rudimentary level, at least some aspects of the 
gas cooling and of the gas transition from hot/warm to 
cold). Questions about the mergings of stellar particles 
were pointed out in Sect. 3.2. Other potentially important 
simplification is the neglection of the growth of metallic- 
ity on the stellar mass-loss. This will be remedied in the 
forthcoming research that will include the chemical en- 
richment. 

We have carried out a simple 2D simulation of a disk 
galaxy starting with no gas. Its purpose was essentially 
illustrative. It should demonstrate the functionality and 
potential applicability of the developed code. Our inten- 
tion was also to show, how the dynamics of classical purely 
stellar simulations could be altered if stellar mass- loss were 
taken into account. 

It goes without saying that the choice of special initial 
conditions - on one hand, zero gas mass in the beginning of 
the simulation and, on the other hand, presumed constant 
star formation rate prior to the simulation - is not coher- 
ent. Nevertheless, such conditions can be approached in 
some astrophysical situations: during strong interactions 
between galaxies, gas is stripped from the outer parts of 
galactic disks (e.g. Moore et al. 1996) while that falling 
towards the center can be largely consumed by a nuclear 
starburst (e.g. Mihos & Hernquist 1996). A quietly star 
forming galaxy can thus be turned into a state without 
much gas. It can however be replenished by the stellar 
mass-loss as is the case of our simple model. Moreover, 
our results as for the rate of mass inflow, enhanced by the 
time-dependent recycling (relative to pure N-body calcu- 
lations), seem not to qualitatively depend on our special 
initial conditions since they point in the same direction 
as the results we obtained for two-component disk models 
(Jungwiert 1998). 

This work will be followed by a second paper dealing 
with the application of the presented code to disk galaxies 
with the emphasis on comparing the instantaneous recy- 
cling approximation with our time-dependent recycling, 
especially as to effects on star formation rate and radial 
transport of matter. 
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